# -*- coding: utf-8 -*-  
'''
lda降维，线性判别分析
    注：仅仅是练手。并没有对大数据做batch处理
    对于小数据量可以做到秒出结果，但大数据量算不出来。。。

Created on 2021年9月16日
@author: luoyi
'''
import numpy as np


#    LDA降维
class LDA:
    '''LDA降维，带标签的降维
        目标：降维后的数据具有一下特征：
            1 类内间距最小。最小化类内方差
            2 类间间距最大。最大化各类均值和
            总结就是保持原有分类
        推导过程：
            假设：K为类别数
                 X为待降维的数据集：[M,N] M为数据维度，N为数据总数
                 Y为降维后的数据：[T,N] K为降维后数据维度
                     Y由X通过线性变换得到：Y = WX 
                     W为变换矩阵：[T, M]
            Y的每个分类均值：
                μ[k,i] = 1/N * ∑(Y∈k)y 
                       = 1/N * ∑(X∈k)W[i]x 
                       = W[i] * 1/N * ∑(X∈k)x
                       = W[i] * μx
                μ[k] = W * μx[k] = Wμx[k]
            类内间距定义为：每个类别的协方差矩阵之和
                cov_y[k] = (Y[k].T - μ[k]) * (Y[k] - μ[k]).T
                         = (WX[k] - Wμx[k])(WX[k] - Wμx[k]).T
                         = W* (X[k] - μx[k])(X[k] - μx[k]).T * W.T
                d_inner = ∑(K) cov_y[K] = ∑(K) W * (X[k] - μx[k])(X[k] - μx[k]).T * W.T
                                        = W * ∑(K)(X[k] - μx[k])(X[k] - μx[k]).T * W.T
                设：Sw = ∑(K)(X[k] - μx[k])(X[k] - μx[k]).T
                      = ∑(K) (n[k]-1) * cov_x[k]
                d_inner = W * Sw * W.T
            类间间距定义为：每个类别均值对总均值的差的加权和。权重为该类别在数据中的占比
                d_outter = ∑(K) nk/N * (μ[k] - μ)(μ[k] - μ).T
                         = ∑(K) nk/N * (W*μx[k] - Wμx)(W*μx[k] - W*μx).T
                         = W * (∑(K)n[k]/N*(μx[k] - μx)(μx[k] - μx).T) * W.T
                设：Sb = ∑(K)n[k]/N * (μx[k] - μx)(μx[k] - μx).T
                d_outter = W * Sb * W.T
            类内间距最小，且类间间距最大，可以表示为：argmax d_outter/d_inner = argmax |(w * Sb * W.T)| / |(w * Sw * W.T)|
            
            矩阵没办法直接算，这里取的是行列式的值（能算就行，别在意这些细节。。。）
                L = argmax(W) |W * Sb * W.T| / |W * Sw * W.T|
                设：W * Sw * W.T = 1
                问题转化为：argmax(W) |W * Sb * W.T|
                                s.t |W * Sw * W.T| = 1
                L(W,λ) = |W * Sb * W.T| - λ(|W * Sw * W.T| - 1)
                ∂L/∂W = 2*Sb*W - 2λ*Sw*W = 0
                      => Sb*W = λ*Sw*W
                      => (Sb*Sw.I) * W = λW
                      => λ为(Sb*Sw.I)的特征值
                         W为其对应的特征向量
                         
            => W即为(Sb*Sw.I)的特征值对应的特征向量，取其最大的前K个特征向量对应的特征值组成的矩阵即为变换矩阵
            
            (Sb*Sw.I)的秩最多是K个，至多有K个特征值，那么降维后至多有K维
            因为 Sb = ∑(K) n[k] * (μ[k] - μ)(μ[k] - μ).T
            ∵ R(AB) <= min(R(A), R(B))，这里(μ[k] - μ)为1维向量
            ∴ 所以对Sb而言每个R((μ[k] - μ)(μ[k] - μ).T) = 1
            ∵ R(A + B) <= R(A) + R(B)
            ∴ 对Sb而言R(Sb) <= ∑(K) ((μ[k] - μ)(μ[k] - μ).T) 
            => R(Sb) <= K
            对Sw.I而言，他都可逆了那他的秩就等于其阶数，Sw的阶数就是原数据维度M
            如果M < K，则可降至M一下任意维度
            如果M > K，则最多只能保留K维
    '''
    def __init__(self, K):
        self._K = K
        pass
    
    #    降维
    def reduce(self, X, labels):
        '''
            @param X: Ndarray(N, M)    需要降维的数据。N：数据长度，M：数据维度
            @param labels: Ndarray(N,) 每个样本所属分类索引
            @return: Y Ndarray(N, K)    降维后的数据
        '''
        #    数据维度
        _ = X.shape[-1]
        #    分类数量，(Sb*Sw.I)特征值（降序），特征值对应的特征向量
        k, _, p = self.w(X, labels)
        #    取前K个特征向量，并规定最大K值=k-1
        K = k-1 if self._K > k-1 else self._K
        w = p[:K+1, :].T
        
        return np.dot(X, w)
    
    def w(self, X, labels):
        d = X.shape[-1]
        N = len(X)                              #    数据总量
        m = np.mean(X, axis=0)                  #    索引数据平均值
        m = np.mat(m)
        
        #    每个分类的索引
        u = np.unique(labels)
        
        #    Sw = ∑(K) (n[k]-1) * cov_x[k]
        Sw = np.mat(np.zeros(shape=(d, d)), dtype=np.float32)
        #    Sb = ∑(K) n[k]/N * (μx[k] - μx)(μx[k] - μx).T
        Sb = np.mat(np.zeros(shape=(d, d)), dtype=np.float32)
        #    计算Sw
        for k in u:
            #    第k个分类的集合
            x = X[np.where(labels == k)[0]]
            
            #    第k个分类的数量
            num_k = len(x)
            
            #    第k个分类的Sw：(n[k]-1) * cov_k
            Sw_k = (num_k-1) * np.cov(x.T)          #    Ndarray(d, d)
            Sw = Sw + Sw_k
            
            #    第k个分类的Sb：n[k]/N * (μx[k] - μx)(μx[k] - μx).T
            m_k = np.mean(x, axis=0)                #    μx[k]
            m_k = np.mat(m_k)
            m_k = m_k - m
            Sb_k = (num_k/N) * np.dot(m_k.T, m_k)   #    Ndarray(d, d)
            Sb = Sb + Sb_k
            pass
        
        #    w = Sb * Sw.I
        S = Sb * np.linalg.inv(Sw)
        
        l, p = np.linalg.eig(S)
        return len(u), l, p
        
    pass
